from numpy import *


def funcion_esfera(x, y, z, R):
    if x**2 + y**2 + z**2 <= R**2:
        return 1
    else:
        return 0

def volumen_esfera(R, N):

    x = random.uniform(-R, R, N)
    y = random.uniform(-R, R, N)
    z = random.uniform(-R, R, N)
    suma = 0
    suma2 = 0

    for i in range(N):
        r = funcion_esfera(x[i], y[i], z[i], R)
        suma += r
        suma2 += r**2

    Vcubo = (2*R)**3
    volumen = Vcubo * suma/N
    error = Vcubo * sqrt((suma2/N - (suma/N)**2)/N)

    return volumen, error

R = 1.66

print("N", "I", "Error", sep="                      ")
print("-"*55)
for n in range(2, 9):
    N = 10**n
    V, error = volumen_esfera(R, N)
    print(f"{N:>12}\t{V:20.12f}\t{error:20.2e}")
print("-"*55)